
###############
## Load data ##
###############

trdmat <- read.csv("data/datasets/main_data.csv")
trdmat$hs10code[nchar(trdmat$hs10code) == 9] <- paste("0", trdmat$hs10code[nchar(trdmat$hs10code) == 9], sep = "")
excl_agg <- read.csv("data/datasets/excl_agg.csv")
tars_covered <- read.csv("data/datasets/tars_covered.csv")

######################################
## Models imports 20122016-20182020 ##
######################################

## top 10
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code, impoth = trdmat$imptop1020122016,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20182020/3, trdmat$impALL20182020/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$imptop10pos20122016, trdmat$imptop10pos20122016, trdmat$imptop10pos20122016, trdmat$imptop10pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_top10 <- newdat2
mod <- lm(log(imp+1) ~ year*market*I(impothpos/10), data = newdat2); 
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1top10 <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2top10 <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## top 25
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20182020/3, trdmat$impALL20182020/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$imptop25pos20122016, trdmat$imptop25pos20122016, trdmat$imptop25pos20122016, trdmat$imptop25pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_top25 <- newdat2
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1top25 <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2top25 <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## ALL
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20182020/3, trdmat$impALL20182020/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALL <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALL <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_top10$hs10code, 1, 6)
mods_top10 <- list(coeftest(mod1top10, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2top10, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_top10[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_top10), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top10)){mod <- mods_top10[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1top10)$df[1:2]),sum(summary(mod2top10)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1_top10_20182020 <- tab

hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_top25 <- list(coeftest(mod1top25, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2top25, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_top25[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_top25), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top25)){mod <- mods_top25[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1top25)$df[1:2]),sum(summary(mod2top25)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1_top25_20182020 <- tab

hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALL <- list(coeftest(mod1ALL, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALL, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALL[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1ALL)$df[1:2]),sum(summary(mod2ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1_ALL_20182020 <- tab

################################################
## Illustrating changes in imports from China ##
################################################

## top 10
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code, impoth = trdmat$imptop1020122016,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20182020/3, trdmat$impALL20182020/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$imptop10pos20122016, trdmat$imptop10pos20122016, trdmat$imptop10pos20122016, trdmat$imptop10pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
mod <- lm(log(imp+1) ~ year*market*I(impothpos/10), data = newdat2)

# with China cfs (fn 24)
qs <- quantile(newdat2$impothpos, c(.10,.5,.9)); names(qs) <- c("","","")
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = 0))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = qs[1]))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = qs[2]))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = qs[3]))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))

## top 25
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20182020/3, trdmat$impALL20182020/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$imptop25pos20122016, trdmat$imptop25pos20122016, trdmat$imptop25pos20122016, trdmat$imptop25pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2)

# with China cfs
qs <- quantile(newdat2$impothpos, c(.10,.5,.9)); names(qs) <- c("","","")
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = 0))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = qs[1]))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = qs[2]))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = qs[3]))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))

## ALL
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20182020/3, trdmat$impALL20182020/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2)

# with China cfs
qs <- quantile(newdat2$impothpos, c(.10,.5,.9)); names(qs) <- c("","","")
preds <- predict(mod, newdata = data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = 0))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))
preds <- predict(mod, newdata = data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = qs[1]))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = qs[2]))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))
preds <- predict(mod, data.frame(year = c("20122016","20182020"), market = c("China"), impothpos = qs[3]))
((exp(preds[2])-1) - (exp(preds[1])-1))/((exp(preds[1])-1))

######################################
## Models imports 20122016-20212023 ##
######################################

## top 10
newdat <- data.frame(year = rep(c("20122016","20212023"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code, impoth = trdmat$imptop1020122016,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20212023/3, trdmat$impALL20212023/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$imptop10pos20122016, trdmat$imptop10pos20122016, trdmat$imptop10pos20122016, trdmat$imptop10pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_top10 <- newdat2
mod <- lm(log(imp+1) ~ year*market*I(impothpos/10), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1top10 <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2top10 <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## top 25
newdat <- data.frame(year = rep(c("20122016","20212023"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20212023/3, trdmat$impALL20212023/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$imptop25pos20122016, trdmat$imptop25pos20122016, trdmat$imptop25pos20122016, trdmat$imptop25pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_top25 <- newdat2
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1top25 <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2top25 <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## ALL
newdat <- data.frame(year = rep(c("20122016","20212023"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20212023/3, trdmat$impALL20212023/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALL <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALL <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_top10$hs10code, 1, 6)
mods_top10 <- list(coeftest(mod1top10, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2top10, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_top10[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_top10), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top10)){mod <- mods_top10[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1top10)$df[1:2]),sum(summary(mod2top10)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1_top10_20212023 <- tab

hs6codes <- substr(newdat2_top25$hs10code, 1, 6)
mods_top25 <- list(coeftest(mod1top25, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2top25, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_top25[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_top25), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top25)){mod <- mods_top25[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1top25)$df[1:2]),sum(summary(mod2top25)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1_top25_20212023 <- tab

hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALL <- list(coeftest(mod1ALL, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALL, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALL[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1ALL)$df[1:2]),sum(summary(mod2ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1_ALL_20212023 <- tab

ptable(cbind(rownames(tab), tab1_top10_20182020, tab1_top10_20212023), "paper/tables/Table1top.tex")
ptable(cbind(rownames(tab), tab1_top25_20182020, tab1_top25_20212023), "paper/tables/Table1mid.tex")
ptable(cbind(rownames(tab), tab1_ALL_20182020, tab1_ALL_20212023), "paper/tables/Table1bot.tex")

########################################
## Trade models allies and non-allies ##
########################################

# summary stats on trade with allies and non-allies
sum(trdmat$impALL20122016allies + trdmat$impALL20182020allies)/sum(trdmat$impALL20122016 + trdmat$impALL20182020)
sum(trdmat$impALL20122016nonallies + trdmat$impALL20182020nonallies)/sum(trdmat$impALL20122016 + trdmat$impALL20182020)

## ALL
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5-trdmat$impchn20122016/5, trdmat$impchn20182020/3, trdmat$impALL20182020/3-trdmat$impchn20182020/3)/1000,
        impallies = c(trdmat$impchn20122016/5, trdmat$impALL20122016allies/5, trdmat$impchn20182020/3, trdmat$impALL20182020allies/3)/1000,
        impnonallies = c(trdmat$impchn20122016/5, trdmat$impALL20122016nonallies/5, trdmat$impchn20182020/3, trdmat$impALL20182020nonallies/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016),
        impothposallies = c(trdmat$impALLpos20122016allies, trdmat$impALLpos20122016allies, trdmat$impALLpos20122016allies, trdmat$impALLpos20122016allies),
        impothposnonallies = c(trdmat$impALLpos20122016nonallies, trdmat$impALLpos20122016nonallies, trdmat$impALLpos20122016nonallies, trdmat$impALLpos20122016nonallies)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(impallies+1) ~ year*market*log(impothposallies+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), impallies = log(newdat2$impallies+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothposallies",
  "year_market","year_impothposallies","market_impothposallies","year_market_impothposallies",
  "impallies","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothposallies_mean",
  "year_market_mean","year_impothposallies_mean","market_impothposallies_mean","year_market_impothposallies_mean",
  "impallies_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALLallies <- lm(impallies ~ year + year_market + year_impothposallies + year_market_impothposallies, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALLallies <- lm(impallies ~ year_impothposallies + year_market_impothposallies + year_sector, data= mm_df); 

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(impnonallies+1) ~ year*market*log(impothposnonallies+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), impnonallies = log(newdat2$impnonallies+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothposnonallies",
  "year_market","year_impothposnonallies","market_impothposnonallies","year_market_impothposnonallies",
  "impnonallies","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothposnonallies_mean",
  "year_market_mean","year_impothposnonallies_mean","market_impothposnonallies_mean","year_market_impothposnonallies_mean",
  "impnonallies_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALLnonallies <- lm(impnonallies ~ year + year_market + year_impothposnonallies + year_market_impothposnonallies, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALLnonallies <- lm(impnonallies ~ year_impothposnonallies + year_market_impothposnonallies + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALLallies <- list(coeftest(mod1ALLallies, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALLallies, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALLallies[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALLallies), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALLallies)){mod <- mods_ALLallies[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Allied alternative markets") 
N <- c(sum(summary(mod1ALLallies)$df[1:2]),sum(summary(mod2ALLallies)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1 <- tab

hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALLnonallies <- list(coeftest(mod1ALLnonallies, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALLnonallies, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALLnonallies[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALLnonallies), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALLnonallies)){mod <- mods_ALLnonallies[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Non-allied alternative markets") 
N <- c(sum(summary(mod1ALLnonallies)$df[1:2]),sum(summary(mod2ALLnonallies)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab2 <- tab

mrg_rw <- rownames(tab); mrg_rw[1] <-  "Post-Trade War$\\cdot$China$\\cdot$Alternative markets"
tab2_ALL_20182020_allies <- tab1
tab2_ALL_20182020_nonallies <- tab2

##########################################
## Trade models friends and non-friends ##
##########################################

# summary stats on trade with friends and non-friends
sum(trdmat$impALL20122016friends + trdmat$impALL20182020friends)/sum(trdmat$impALL20122016 + trdmat$impALL20182020)
sum(trdmat$impALL20122016nonfriends + trdmat$impALL20182020nonfriends)/sum(trdmat$impALL20122016 + trdmat$impALL20182020)

## ALL
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20182020/3, trdmat$impALL20182020/3)/1000,
        impfriends = c(trdmat$impchn20122016/5, trdmat$impALL20122016friends/5, trdmat$impchn20182020/3, trdmat$impALL20182020friends/3)/1000,
        impnonfriends = c(trdmat$impchn20122016/5, trdmat$impALL20122016nonfriends/5, trdmat$impchn20182020/3, trdmat$impALL20182020nonfriends/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016),
        impothposfriends = c(trdmat$impALLpos20122016friends, trdmat$impALLpos20122016friends, trdmat$impALLpos20122016friends, trdmat$impALLpos20122016friends),
        impothposnonfriends = c(trdmat$impALLpos20122016nonfriends, trdmat$impALLpos20122016nonfriends, trdmat$impALLpos20122016nonfriends, trdmat$impALLpos20122016nonfriends)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(impfriends+1) ~ year*market*log(impothposfriends+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), impfriends = log(newdat2$impfriends+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothposfriends",
  "year_market","year_impothposfriends","market_impothposfriends","year_market_impothposfriends",
  "impfriends","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothposfriends_mean",
  "year_market_mean","year_impothposfriends_mean","market_impothposfriends_mean","year_market_impothposfriends_mean",
  "impfriends_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALLfriends <- lm(impfriends ~ year + year_market + year_impothposfriends + year_market_impothposfriends, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALLfriends <- lm(impfriends ~ year_impothposfriends + year_market_impothposfriends + year_sector, data= mm_df); 

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(impnonfriends+1) ~ year*market*log(impothposnonfriends+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), impnonfriends = log(newdat2$impnonfriends+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothposnonfriends",
  "year_market","year_impothposnonfriends","market_impothposnonfriends","year_market_impothposnonfriends",
  "impnonfriends","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothposnonfriends_mean",
  "year_market_mean","year_impothposnonfriends_mean","market_impothposnonfriends_mean","year_market_impothposnonfriends_mean",
  "impnonfriends_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALLnonfriends <- lm(impnonfriends ~ year + year_market + year_impothposnonfriends + year_market_impothposnonfriends, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALLnonfriends <- lm(impnonfriends ~ year_impothposnonfriends + year_market_impothposnonfriends + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALLfriends <- list(coeftest(mod1ALLfriends, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALLfriends, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALLfriends[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALLfriends), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALLfriends)){mod <- mods_ALLfriends[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Friendly alternative markets") 
N <- c(sum(summary(mod1ALLfriends)$df[1:2]),sum(summary(mod2ALLfriends)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1 <- tab

hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALLnonfriends <- list(coeftest(mod1ALLnonfriends, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALLnonfriends, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALLnonfriends[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALLnonfriends), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALLnonfriends)){mod <- mods_ALLnonfriends[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Non-friendly alternative markets") 
N <- c(sum(summary(mod1ALLnonfriends)$df[1:2]),sum(summary(mod2ALLnonfriends)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab2 <- tab

mrg_rw <- rownames(tab); mrg_rw[1] <-  "Post-Trade War$\\cdot$China$\\cdot$Alternative markets"
tab2_ALL_20182020_friends <- tab1
tab2_ALL_20182020_nonfriends <- tab2

#################################################
## Trade models allies and non-allies 20212023 ##
#################################################

## ALL
newdat <- data.frame(year = rep(c("20122016","20212023"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20212023/3, trdmat$impALL20212023/3)/1000,
        impallies = c(trdmat$impchn20122016/5, trdmat$impALL20122016allies/5, trdmat$impchn20212023/3, trdmat$impALL20212023allies/3)/1000,
        impnonallies = c(trdmat$impchn20122016/5, trdmat$impALL20122016nonallies/5, trdmat$impchn20212023/3, trdmat$impALL20212023nonallies/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016),
        impothposallies = c(trdmat$impALLpos20122016allies, trdmat$impALLpos20122016allies, trdmat$impALLpos20122016allies, trdmat$impALLpos20122016allies),
        impothposnonallies = c(trdmat$impALLpos20122016nonallies, trdmat$impALLpos20122016nonallies, trdmat$impALLpos20122016nonallies, trdmat$impALLpos20122016nonallies)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(impallies+1) ~ year*market*log(impothposallies+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), impallies = log(newdat2$impallies+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothposallies",
  "year_market","year_impothposallies","market_impothposallies","year_market_impothposallies",
  "impallies","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothposallies_mean",
  "year_market_mean","year_impothposallies_mean","market_impothposallies_mean","year_market_impothposallies_mean",
  "impallies_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALLallies <- lm(impallies ~ year + year_market + year_impothposallies + year_market_impothposallies, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALLallies <- lm(impallies ~ year_impothposallies + year_market_impothposallies + year_sector, data= mm_df); 

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(impnonallies+1) ~ year*market*log(impothposnonallies+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), impnonallies = log(newdat2$impnonallies+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothposnonallies",
  "year_market","year_impothposnonallies","market_impothposnonallies","year_market_impothposnonallies",
  "impnonallies","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothposnonallies_mean",
  "year_market_mean","year_impothposnonallies_mean","market_impothposnonallies_mean","year_market_impothposnonallies_mean",
  "impnonallies_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALLnonallies <- lm(impnonallies ~ year + year_market + year_impothposnonallies + year_market_impothposnonallies, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALLnonallies <- lm(impnonallies ~ year_impothposnonallies + year_market_impothposnonallies + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALLallies <- list(coeftest(mod1ALLallies, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALLallies, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALLallies[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALLallies), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALLallies)){mod <- mods_ALLallies[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Allied alternative markets") 
N <- c(sum(summary(mod1ALLallies)$df[1:2]),sum(summary(mod2ALLallies)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1 <- tab

hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALLnonallies <- list(coeftest(mod1ALLnonallies, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALLnonallies, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALLnonallies[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALLnonallies), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALLnonallies)){mod <- mods_ALLnonallies[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Non-allied alternative markets") 
N <- c(sum(summary(mod1ALLnonallies)$df[1:2]),sum(summary(mod2ALLnonallies)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab2 <- tab

mrg_rw <- rownames(tab); mrg_rw[1] <-  "Post-Trade War$\\cdot$China$\\cdot$Alternative markets"
tab2_ALL_20212023_allies <- tab1
tab2_ALL_20212023_nonallies <- tab2

ptable(cbind(mrg_rw, tab2_ALL_20182020_allies, tab2_ALL_20212023_allies), "paper/tables/Table2toptop.tex")
ptable(cbind(mrg_rw, tab2_ALL_20182020_nonallies, tab2_ALL_20212023_nonallies), "paper/tables/Table2topbot.tex")

###################################################
## Trade models friends and non-friends 20212023 ##
###################################################

## ALL
newdat <- data.frame(year = rep(c("20122016","20212023"), each = 2*nrow(trdmat)), 
        hs10code = trdmat$hs10code,
        imp = c(trdmat$impchn20122016/5, trdmat$impALL20122016/5, trdmat$impchn20212023/3, trdmat$impALL20212023/3)/1000,
        impfriends = c(trdmat$impchn20122016/5, trdmat$impALL20122016friends/5, trdmat$impchn20212023/3, trdmat$impALL20212023friends/3)/1000,
        impnonfriends = c(trdmat$impchn20122016/5, trdmat$impALL20122016nonfriends/5, trdmat$impchn20212023/3, trdmat$impALL20212023nonfriends/3)/1000,
        market = c(rep("China", nrow(trdmat)), rep("Other", nrow(trdmat)), rep("China", nrow(trdmat)), rep("Other", nrow(trdmat))),
        impothpos = c(trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016, trdmat$impALLpos20122016),
        impothposfriends = c(trdmat$impALLpos20122016friends, trdmat$impALLpos20122016friends, trdmat$impALLpos20122016friends, trdmat$impALLpos20122016friends),
        impothposnonfriends = c(trdmat$impALLpos20122016nonfriends, trdmat$impALLpos20122016nonfriends, trdmat$impALLpos20122016nonfriends, trdmat$impALLpos20122016nonfriends)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(impfriends+1) ~ year*market*log(impothposfriends+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), impfriends = log(newdat2$impfriends+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothposfriends",
  "year_market","year_impothposfriends","market_impothposfriends","year_market_impothposfriends",
  "impfriends","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothposfriends_mean",
  "year_market_mean","year_impothposfriends_mean","market_impothposfriends_mean","year_market_impothposfriends_mean",
  "impfriends_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALLfriends <- lm(impfriends ~ year + year_market + year_impothposfriends + year_market_impothposfriends, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALLfriends <- lm(impfriends ~ year_impothposfriends + year_market_impothposfriends + year_sector, data= mm_df); 

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(impnonfriends+1) ~ year*market*log(impothposnonfriends+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), impnonfriends = log(newdat2$impnonfriends+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothposnonfriends",
  "year_market","year_impothposnonfriends","market_impothposnonfriends","year_market_impothposnonfriends",
  "impnonfriends","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs10code_market <- paste(mm_df$hs10code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs10code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code_market","year_mean","market_mean","impothposnonfriends_mean",
  "year_market_mean","year_impothposnonfriends_mean","market_impothposnonfriends_mean","year_market_impothposnonfriends_mean",
  "impnonfriends_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALLnonfriends <- lm(impnonfriends ~ year + year_market + year_impothposnonfriends + year_market_impothposnonfriends, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALLnonfriends <- lm(impnonfriends ~ year_impothposnonfriends + year_market_impothposnonfriends + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALLfriends <- list(coeftest(mod1ALLfriends, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALLfriends, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALLfriends[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALLfriends), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALLfriends)){mod <- mods_ALLfriends[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Friendly alternative markets") 
N <- c(sum(summary(mod1ALLfriends)$df[1:2]),sum(summary(mod2ALLfriends)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1 <- tab

hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALLnonfriends <- list(coeftest(mod1ALLnonfriends, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALLnonfriends, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALLnonfriends[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALLnonfriends), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALLnonfriends)){mod <- mods_ALLnonfriends[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Non-friendly alternative markets") 
N <- c(sum(summary(mod1ALLnonfriends)$df[1:2]),sum(summary(mod2ALLnonfriends)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab2 <- tab

mrg_rw <- rownames(tab); mrg_rw[1] <-  "Post-Trade War$\\cdot$China$\\cdot$Alternative markets"
tab2_ALL_20212023_friends <- tab1
tab2_ALL_20212023_nonfriends <- tab2

ptable(cbind(mrg_rw, tab2_ALL_20182020_friends, tab2_ALL_20212023_friends), "paper/tables/Table2bottop.tex")
ptable(cbind(mrg_rw, tab2_ALL_20182020_nonfriends, tab2_ALL_20212023_nonfriends), "paper/tables/Table2botbot.tex")

############################### 
## Models exclusion requests ##
###############################

## top 10
newdat <- data.frame(hs10code = trdmat$hs10code,
        impchn = c(trdmat$impchn20122016/5)/1000,
        imptop10 = c(trdmat$imptop1020122016/5)/1000,
        imptop10pos = c(trdmat$imptop10pos20122016),
        impwornonchnpos = trdmat$impALL20122016 > 0, 
        impwor = c(trdmat$impALL20122016/5/1000), 
        TWtariff = trdmat$TWtariff)

# merge in exclusion requestions variable
newdat <- join(newdat, excl_agg, by = "hs10code", type = "left")
newdat$numexclreq[is.na(newdat$numexclreq)] <- 0

newdat$covered <- 0
newdat$covered[substr(newdat$hs10code,1,8) %in% tars_covered$hs10code] <- 1

ana_dat <- na.omit(newdat[,c("numexclreq","imptop10pos","impwor","impchn","covered","hs10code")])
postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
ana_dat <- ana_dat[ana_dat$hs10code %in% postrade,]

mod1a_top10 <- lm(log(numexclreq+1) ~ I(imptop10pos/10) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat);
mod1b_top10 <- lm(log(numexclreq+1) ~ I(imptop10pos/10) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

cf0 <- cf1 <- ana_dat; cf0$imptop10pos <- 0; cf1$imptop10pos <- 9
mean(predict(mod1a_top10, cf1)) - mean(predict(mod1a_top10, cf0))

mod2a_top10 <- lm(I(numexclreq>0) ~ I(imptop10pos/10) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod2b_top10 <- lm(I(numexclreq>0) ~ I(imptop10pos/10) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

cf0 <- cf1 <- ana_dat; cf0$imptop10pos <- 0; cf1$imptop10pos <- 9
mean(predict(mod2a_top10, cf1)) - mean(predict(mod2a_top10, cf0))

mods_top10 <- list(coeftest(mod1a_top10, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod1b_top10, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2a_top10, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2b_top10, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4]) 

vars <- c(rownames(mods_top10[[1]])[c(2)])
tab <- matrix(data = NA, ncol = length(mods_top10), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top10)){mod <- mods_top10[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Alternative markets") 
N <- c(sum(summary(mod1a_top10)$df[1:2]),sum(summary(mod1b_top10)$df[1:2]),sum(summary(mod2a_top10)$df[1:2]),sum(summary(mod2b_top10)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))

ptable(cbind(rownames(tab), tab)[,1:3], "paper/tables/Table3top.tex")
ptable(cbind(rownames(tab), tab)[,c(1,4:5)], "paper/tables/TableA3top.tex")

## top 25
newdat <- data.frame(hs10code = trdmat$hs10code,
        impchn = c(trdmat$impchn20122016/5)/1000,
        imptop25 = c(trdmat$imptop2520122016/5)/1000,
        imptop25pos = c(trdmat$imptop25pos20122016),
        impwornonchnpos = trdmat$impALL20122016 > 0,
        impwor = c(trdmat$impALL20122016/5/1000),
        TWtariff = trdmat$TWtariff)

# merge in exclusion requestions variable
newdat <- join(newdat, excl_agg, by = "hs10code", type = "left")
newdat$numexclreq[is.na(newdat$numexclreq)] <- 0

newdat$covered <- 0
newdat$covered[substr(newdat$hs10code,1,8) %in% tars_covered$hs10code] <- 1

ana_dat <- na.omit(newdat[,c("numexclreq","imptop25pos","impwor","impchn","covered","hs10code")])
postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
ana_dat <- ana_dat[ana_dat$hs10code %in% postrade,]

mod1a_top25 <- lm(log(numexclreq+1) ~ log(imptop25pos+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod1b_top25 <- lm(log(numexclreq+1) ~ log(imptop25pos+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mod2a_top25 <- lm(I(numexclreq>0) ~ log(imptop25pos+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod2b_top25 <- lm(I(numexclreq>0) ~ log(imptop25pos+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mods_top25 <- list(coeftest(mod1a_top25, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod1b_top25, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2a_top25, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2b_top25, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4]) 

vars <- c(rownames(mods_top25[[1]])[c(2)])
tab <- matrix(data = NA, ncol = length(mods_top25), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top25)){mod <- mods_top25[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Alternative markets") 
N <- c(sum(summary(mod1a_top25)$df[1:2]),sum(summary(mod1b_top25)$df[1:2]),sum(summary(mod2a_top25)$df[1:2]),sum(summary(mod2b_top25)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))

ptable(cbind(rownames(tab), tab)[,1:3], "paper/tables/Table3mid.tex")
ptable(cbind(rownames(tab), tab)[,c(1,4:5)], "paper/tables/TableA3mid.tex")

## ALL
newdat <- data.frame(hs10code = trdmat$hs10code,
        impchn = c(trdmat$impchn20122016/5)/1000,
        impALL = c(trdmat$impALL20122016/5)/1000,
        impALLpos = c(trdmat$impALLpos20122016),
        impwornonchnpos = trdmat$impALL20122016 > 0,
        impwor = c(trdmat$impALL20122016/5/1000),
        TWtariff = trdmat$TWtariff)

# merge in exclusion requestions variable
newdat <- join(newdat, excl_agg, by = "hs10code", type = "left")
newdat$numexclreq[is.na(newdat$numexclreq)] <- 0

newdat$covered <- 0
newdat$covered[substr(newdat$hs10code,1,8) %in% tars_covered$hs10code] <- 1

ana_dat <- na.omit(newdat[,c("numexclreq","impALLpos","impwor","impchn","covered","hs10code","TWtariff")])
postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
ana_dat <- ana_dat[ana_dat$hs10code %in% postrade,]

mod1a_ALL <- lm(log(numexclreq+1) ~ log(impALLpos+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod1b_ALL <- lm(log(numexclreq+1) ~ log(impALLpos+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mod2a_ALL <- lm(I(numexclreq>0) ~ log(impALLpos+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod2b_ALL <- lm(I(numexclreq>0) ~ log(impALLpos+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mods_ALL <- list(coeftest(mod1a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod1b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4]) 

vars <- c(rownames(mods_ALL[[1]])[c(2)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Alternative markets") 
N <- c(sum(summary(mod1a_ALL)$df[1:2]),sum(summary(mod1b_ALL)$df[1:2]),sum(summary(mod2a_ALL)$df[1:2]),sum(summary(mod2b_ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))

mrg_rw <- rownames(tab); mrg_rw[1] <- "Alternative markets"
ptable(cbind(mrg_rw, tab)[,1:3], "paper/tables/Table3bot.tex")
ptable(cbind(rownames(tab), tab)[,c(1,4:5)], "paper/tables/TableA3bot.tex")

####################################################
## Models exclusion requests allies and nonallies ##
####################################################

## ALL
newdat <- data.frame(hs10code = trdmat$hs10code,
        impchn = c(trdmat$impchn20122016/5)/1000,
        impALL = c(trdmat$impALL20122016/5)/1000,
        impALLposallies = c(trdmat$impALLpos20122016allies),
        impALLposnonallies = c(trdmat$impALLpos20122016nonallies),
        impwornonchnpos = trdmat$impALL20122016 > 0,
        impwor = c(trdmat$impALL20122016/5/1000),
        TWtariff = trdmat$TWtariff)

# merge in exclusion requestions variable
newdat <- join(newdat, excl_agg, by = "hs10code", type = "left")
newdat$numexclreq[is.na(newdat$numexclreq)] <- 0

newdat$covered <- 0
newdat$covered[substr(newdat$hs10code,1,8) %in% tars_covered$hs10code] <- 1

ana_dat <- na.omit(newdat[,c("numexclreq","impALLposallies","impwor","impchn","covered","hs10code","TWtariff")])
postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
ana_dat <- ana_dat[ana_dat$hs10code %in% postrade,]

mod1a_ALL <- lm(log(numexclreq+1) ~ log(impALLposallies+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod1b_ALL <- lm(log(numexclreq+1) ~ log(impALLposallies+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mod2a_ALL <- lm(I(numexclreq>0) ~ log(impALLposallies+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod2b_ALL <- lm(I(numexclreq>0) ~ log(impALLposallies+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mods_ALL <- list(coeftest(mod1a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod1b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4]) 

vars <- c(rownames(mods_ALL[[1]])[c(2)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Alternative markets") 
N <- c(sum(summary(mod1a_ALL)$df[1:2]),sum(summary(mod1b_ALL)$df[1:2]),sum(summary(mod2a_ALL)$df[1:2]),sum(summary(mod2b_ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1 <- tab[,1:2]
ptable(cbind(rownames(tab1), tab1), "paper/tables/Table4toptop.tex")

ana_dat <- na.omit(newdat[,c("numexclreq","impALLposnonallies","impwor","impchn","covered","hs10code","TWtariff")])
postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
ana_dat <- ana_dat[ana_dat$hs10code %in% postrade,]

mod1a_ALL <- lm(log(numexclreq+1) ~ log(impALLposnonallies+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod1b_ALL <- lm(log(numexclreq+1) ~ log(impALLposnonallies+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mod2a_ALL <- lm(I(numexclreq>0) ~ log(impALLposnonallies+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat);
mod2b_ALL <- lm(I(numexclreq>0) ~ log(impALLposnonallies+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mods_ALL <- list(coeftest(mod1a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod1b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4]) 

vars <- c(rownames(mods_ALL[[1]])[c(2)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Alternative markets") 
N <- c(sum(summary(mod1a_ALL)$df[1:2]),sum(summary(mod1b_ALL)$df[1:2]),sum(summary(mod2a_ALL)$df[1:2]),sum(summary(mod2b_ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab2 <- tab[,1:2]
ptable(cbind(rownames(tab2), tab2), "paper/tables/Table4topbot.tex")

######################################################
## Models exclusion requests friends and nonfriends ##
######################################################

## ALL
newdat <- data.frame(hs10code = trdmat$hs10code,
        impchn = c(trdmat$impchn20122016/5)/1000,
        impALL = c(trdmat$impALL20122016/5)/1000,
        impALLposfriends = c(trdmat$impALLpos20122016friends),
        impALLposnonfriends = c(trdmat$impALLpos20122016nonfriends),
        impwornonchnpos = trdmat$impALL20122016 > 0,
        impwor = c(trdmat$impALL20122016/5/1000),
        TWtariff = trdmat$TWtariff)

# merge in exclusion requestions variable
newdat <- join(newdat, excl_agg, by = "hs10code", type = "left")
newdat$numexclreq[is.na(newdat$numexclreq)] <- 0

newdat$covered <- 0
newdat$covered[substr(newdat$hs10code,1,8) %in% tars_covered$hs10code] <- 1

ana_dat <- na.omit(newdat[,c("numexclreq","impALLposfriends","impwor","impchn","covered","hs10code","TWtariff")])
postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
ana_dat <- ana_dat[ana_dat$hs10code %in% postrade,]

mod1a_ALL <- lm(log(numexclreq+1) ~ log(impALLposfriends+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod1b_ALL <- lm(log(numexclreq+1) ~ log(impALLposfriends+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mod2a_ALL <- lm(I(numexclreq>0) ~ log(impALLposfriends+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat);
mod2b_ALL <- lm(I(numexclreq>0) ~ log(impALLposfriends+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mods_ALL <- list(coeftest(mod1a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod1b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4]) 

vars <- c(rownames(mods_ALL[[1]])[c(2)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Alternative markets") 
N <- c(sum(summary(mod1a_ALL)$df[1:2]),sum(summary(mod1b_ALL)$df[1:2]),sum(summary(mod2a_ALL)$df[1:2]),sum(summary(mod2b_ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1 <- tab[,1:2]

ptable(cbind(rownames(tab1), tab1), "paper/tables/Table4bottop.tex")

ana_dat <- na.omit(newdat[,c("numexclreq","impALLposnonfriends","impwor","impchn","covered","hs10code","TWtariff")])
postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
ana_dat <- ana_dat[ana_dat$hs10code %in% postrade,]

mod1a_ALL <- lm(log(numexclreq+1) ~ log(impALLposnonfriends+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat); 
mod1b_ALL <- lm(log(numexclreq+1) ~ log(impALLposnonfriends+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mod2a_ALL <- lm(I(numexclreq>0) ~ log(impALLposnonfriends+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,2), data = ana_dat);
mod2b_ALL <- lm(I(numexclreq>0) ~ log(impALLposnonfriends+1) + log10(impwor+1) + log10(impchn+1) + covered + substr(hs10code,1,4), data = ana_dat); 

mods_ALL <- list(coeftest(mod1a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod1b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2a_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4], 
  coeftest(mod2b_ALL, vcov = vcovCL, cluster = substr(ana_dat$hs10code,1,6))[1:3,1:4]) 

vars <- c(rownames(mods_ALL[[1]])[c(2)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Alternative markets") 
N <- c(sum(summary(mod1a_ALL)$df[1:2]),sum(summary(mod1b_ALL)$df[1:2]),sum(summary(mod2a_ALL)$df[1:2]),sum(summary(mod2b_ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab2 <- tab[,1:2]
ptable(cbind(rownames(tab2), tab2), "paper/tables/Table4botbot.tex")

############
## Figure ## 
############

## examples of industries with lots of hard to replace China products.
dire <- as.numeric(trdmat$impchn20122016/trdmat$impwor20122016 > .5 & trdmat$imptop10pos20122016 <= 5)
agg <- aggregate(dire ~ substr(trdmat$hs10code, 1, 2), FUN = mean, na.rm = TRUE)
agg2 <- aggregate(rep(1,nrow(trdmat)) ~ substr(trdmat$hs10code, 1, 2), FUN = sum, na.rm = TRUE)
agg <- join(agg, agg2, type = "left")
names(agg) <- c("hs2code","dire","count")

agg <- agg[order(agg$dire, decreasing = TRUE),]
agg[agg$count > 100,][1:10,]

dire <- as.numeric(trdmat$impchn20122016/trdmat$impwor20122016 > .5 & trdmat$impALLpos20122016 <= 10)
agg <- aggregate(dire ~ substr(trdmat$hs10code, 1, 2), FUN = mean, na.rm = TRUE)
agg2 <- aggregate(rep(1,nrow(trdmat)) ~ substr(trdmat$hs10code, 1, 2), FUN = sum, na.rm = TRUE)
agg <- join(agg, agg2, type = "left")
names(agg) <- c("hs2code","dire","count")

agg <- agg[order(agg$dire, decreasing = TRUE),]
agg[agg$count > 100,][1:10,]

##
pdf("paper/figs/Figure2.pdf", height = 6, width = 3.5)
par(mfrow = c(3,1), mar=c(0,1.5,1.25,.5), oma = c(.5, .2, 0, 0), lwd = .5) # bottom, left, top, right

barplot(table(trdmat$imptop10pos20122016[trdmat$impchn20122016 > 0]), axes = F, main = "",
  names.arg = "", ylim = c(-500,7000))
text(x = 1.2*seq(0,10,1) + .7, y = -300, labels = seq(0,10,1), cex = .6)
axis(2, cex.axis = .4, tck = -.015, padj = 3.6, labels = seq(0,6000,1000), at = seq(0,6000,1000))
title(main="# of non-China top-10 markets producing product", cex.main = .6, line = -.5)
mtext("Number of products", 2, line=.85, cex = .5)

barplot(table(trdmat$imptop25pos20122016[trdmat$impchn20122016 > 0]), axes = F, main = "",
  names.arg = "", ylim = c(-80,950))
text(x = 1.2*seq(0,25,1) + .7, y = -47, labels = seq(0,25,1), cex = .5)
axis(2, cex.axis = .4, tck = -.015, padj = 3.6)
title(main="# of non-China top-25 markets producing product", cex.main = .6, line = -.5)
mtext("Number of products", 2, line=.85, cex = .5)

barplot(table(trdmat$impALLpos20122016[trdmat$impchn20122016 > 0]), axes = F, main = "",
  names.arg = "", ylim = c(-80,685), xlim = c(0,200), border = "gray")

dens1 <- dens1orig <- density(trdmat$impALLpos20122016[trdmat$impchn20122016 > 0], n = 460)
dens1$y <- nrow(trdmat)/sum(dens1$y) * dens1$y
dens1$x <- dens1$x*1.2 + .7
mycol1 <- rgb(r = 0, g = 0, b = 0, alpha = .8)
lines(dens1, lwd = 1.2, col = mycol1)

dens2 <- density(trdmat$impALLpos20122016allies[trdmat$impchn20122016 > 0], n = 460)
dens2$y <- nrow(trdmat)/sum(dens1orig$y) * dens2$y
dens2$x <- dens2$x*1.2 + .7
mycol2 <- rgb(r = .66, g = .66, b = .66, alpha = .8)
lines(dens2, lwd = 1.2, col = mycol2)

dens3 <- density(trdmat$impALLpos20122016friends[trdmat$impchn20122016 > 0], n = 460)
dens3$y <- nrow(trdmat)/sum(dens1orig$y) * dens3$y
dens3$x <- dens3$x*1.2 + .7
mycol3 <- rgb(r = .83, g = .83, b = .83, alpha = .8)
lines(dens3, lwd = 1.2, col = mycol3)

text(x = 1.2*seq(0,160,10) + .7, y = -47, labels = seq(0,160,10), cex = .5)
axis(2, cex.axis = .4, tck = -.015, padj = 3.6, labels = seq(0,600,100), at = seq(0,600,100))
title(main="# of non-China markets producing product", cex.main = .6, line = -.5)
mtext("Number of products", 2, line=.85, cex = .5)

legend(x = 125, y = 550, legend = c("All partners","Allied partners","Friendly partners"),
  lwd = 2, col = c(mycol1,mycol2,mycol3), cex = .5, title = "Density curves")

dev.off()

# Figure 1
trd <- read.csv("data/datasets/us china trade.csv", skip = 1)

trd$imp_chn_per <- 100*trd$imp_chn/trd$imp_wor
trd$exp_chn_per <- 100*trd$exp_chn/trd$exp_wor
trd$fdi_chn_per <- 100*trd$fdi_chn/trd$fdi_wor
trd$dia_chn_per <- 100*trd$dia_chn/trd$dia_wor
trd$gdpdef <- trd$gdpdef/100

trd <- trd[trd$year > 1985,]

pdf("paper/figs/Figure1.pdf", height = 7, width = 5)
par(mfrow = c(4,1), mar=c(1,2.5,1.2,1.5), oma = c(.5, .2, 0, 0), lwd = .5) # bottom, left, top, right

plot((imp_chn/1000)/gdpdef ~ year, data = trd, ylim = c(0,550), type = "line", lwd = 2, axes = F)
lines((exp_chn/1000)/gdpdef ~ year, data = trd, lty = "dashed", lwd = 2)
axis(1, cex.axis = .7, tck = -.02, padj = -1.5,
  at = seq(1985,2025,by=5), labels = c(seq(1985,2020,by=5),""))
axis(2, cex.axis = .7, tck = -.015, padj = 1.5)
title(main="Goods trade with China", cex.main = .7, line = -.5)
mtext("Billions of (2012) US dollars", 2, line=1.5, cex = .5)
legend(x = 1986, y = 500,lty = c("solid","dashed"), legend = c("US imports","US exports"),
  lwd = 1, cex = .64)

plot(imp_chn_per ~ year, data = trd, ylim = c(0,25), type = "line", lwd = 2, axes = F, 
  main = "", xlab = "", ylab = "")
lines(exp_chn_per ~ year, data = trd, lty = "dashed", lwd = 2)
axis(1, cex.axis = .7, tck = -.015, padj = -1.5,
  at = seq(1985,2025,by=5), labels = c(seq(1985,2020,by=5),""))
axis(2, cex.axis = .7, tck = -.02, padj = 1.5)
title(main="Goods trade with China (% of US total)", cex.main = .7, line = -.5)
mtext("Percentage", 2, line=1.5, cex = .5)
legend(x = 1986, y = 25,lty = c("solid","dashed"), legend = c("US imports","US exports"),
  lwd = 1, cex = .64)

plot((dia_chn/1000)/gdpdef ~ year, data = trd, ylim = c(0,150), type = "line", lwd = 2, axes = F)
lines((fdi_chn/1000)/gdpdef ~ year, data = trd, lty = "dashed", lwd = 2)
axis(1, cex.axis = .7, tck = -.015, padj = -1.5,
  at = seq(1985,2025,by=5), labels = c(seq(1985,2020,by=5),""))
axis(2, cex.axis = .7, tck = -.02, padj = 1.5)
title(main="FDI stocks with China", cex.main = .7, line = -.5)
mtext("Billions of (2012) US dollars", 2, line=1.5, cex = .5)
legend(x = 1986, y = 150,lty = c("solid","dashed"), legend = c("US outward","US inward"),
  lwd = 1, cex = .64)

plot(dia_chn_per ~ year, data = trd, ylim = c(0,5), type = "line", lwd = 2, axes = F)
lines(fdi_chn_per ~ year, data = trd, lty = "dashed", lwd = 2)
axis(1, cex.axis = .7, tck = -.015, padj = -1.5,
  at = seq(1985,2025,by=5), labels = c(seq(1985,2020,by=5),""))
axis(2, cex.axis = .7, tck = -.02, padj = 1.5)
title(main="FDI stocks with China (% of US total)", cex.main = .7, line = -.5)
# mtext("Year", 1, line=1, cex = .5)
mtext("Percentage", 2, line=1.5, cex = .5)
legend(x = 1986, y = 5, lty = c("solid","dashed"), legend = c("US outward","US inward"),
  lwd = 1, cex = .64)

dev.off()

##############################
## Models imports TOT TRADE ##
##############################

### 20122016 to 20182020
## top 10
newdat <- data.frame(year = rep(c("20122016","20182020"), each = nrow(trdmat)), 
        hs10code = trdmat$hs10code, impoth = trdmat$imptop1020122016,
        imp = c(trdmat$impchn20122016/5, trdmat$impchn20182020/3)/1000
)

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_top10 <- newdat2
mod <- lm(log(imp+1) ~ year*log(impoth+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","impoth","year_impoth","imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:5]], 2, as.numeric)
agg <- aggregate(aggvars ~ hs10code, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code","year_mean","impoth_mean",
  "year_impoth_mean","imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:5] <- mm_df[,2:5] - mm_df[,7:10] 

mod1top10 <- lm(imp ~ year + year_impoth, data = mm_df)
mm_df$year_sector <- paste(newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2top10 <- lm(imp ~ year + year_impoth + year_sector, data= mm_df); 

## top 25
newdat <- data.frame(year = rep(c("20122016","20182020"), each = nrow(trdmat)), 
        hs10code = trdmat$hs10code, impoth = trdmat$imptop2520122016,
        imp = c(trdmat$impchn20122016/5, trdmat$impchn20182020/3)/1000
)

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_top25 <- newdat2
mod <- lm(log(imp+1) ~ year*log(impoth+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","impoth","year_impoth","imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:5]], 2, as.numeric)
agg <- aggregate(aggvars ~ hs10code, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code","year_mean","impoth_mean",
  "year_impoth_mean","imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:5] <- mm_df[,2:5] - mm_df[,7:10] 

mod1top25 <- lm(imp ~ year + year_impoth, data = mm_df)
mm_df$year_sector <- paste(newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2top25 <- lm(imp ~ year + year_impoth + year_sector, data= mm_df); 

## ALL
newdat <- data.frame(year = rep(c("20122016","20182020"), each = nrow(trdmat)), 
        hs10code = trdmat$hs10code, impoth = trdmat$impALL20122016,
        imp = c(trdmat$impchn20122016/5, trdmat$impchn20182020/3)/1000
)

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(imp+1) ~ year*log(impoth+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","impoth","year_impoth","imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:5]], 2, as.numeric)
agg <- aggregate(aggvars ~ hs10code, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code","year_mean","impoth_mean",
  "year_impoth_mean","imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:5] <- mm_df[,2:5] - mm_df[,7:10] 

mod1ALL <- lm(imp ~ year + year_impoth, data = mm_df)
mm_df$year_sector <- paste(newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALL <- lm(imp ~ year + year_impoth + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_top10$hs10code, 1, 6)
mods_top10 <- list(coeftest(mod1top10, vcov = vcovCL, cluster = hs6codes)[1:3,1:4], 
  coeftest(mod2top10, vcov = vcovCL, cluster = hs6codes)[1:3,1:4])

vars <- c(rownames(mods_top10[[1]])[c(3)])
tab <- matrix(data = NA, ncol = length(mods_top10), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top10)){mod <- mods_top10[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$Alternative markets total trade") 
N <- c(sum(summary(mod1top10)$df[1:2]),sum(summary(mod2top10)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1820_top10 <- tab

hs6codes <- substr(newdat2_top25$hs10code, 1, 6)
mods_top25 <- list(coeftest(mod1top25, vcov = vcovCL, cluster = hs6codes)[1:3,1:4], 
  coeftest(mod2top25, vcov = vcovCL, cluster = hs6codes)[1:3,1:4])

vars <- c(rownames(mods_top25[[1]])[c(3)])
tab <- matrix(data = NA, ncol = length(mods_top25), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top25)){mod <- mods_top25[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$Alternative markets total trade") 
N <- c(sum(summary(mod1top25)$df[1:2]),sum(summary(mod2top25)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1820_top25 <- tab

hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALL <- list(coeftest(mod1ALL, vcov = vcovCL, cluster = hs6codes)[1:3,1:4], 
  coeftest(mod2ALL, vcov = vcovCL, cluster = hs6codes)[1:3,1:4])

vars <- c(rownames(mods_ALL[[1]])[c(3)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$Alternative markets total trade") 
N <- c(sum(summary(mod1ALL)$df[1:2]),sum(summary(mod2ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab1820_ALL <- tab

### 20122016 to 20212023
## top 10
newdat <- data.frame(year = rep(c("20122016","20212023"), each = nrow(trdmat)), 
        hs10code = trdmat$hs10code, impoth = trdmat$imptop1020122016,
        imp = c(trdmat$impchn20122016/5, trdmat$impchn20212023/3)/1000
)

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_top10 <- newdat2
mod <- lm(log(imp+1) ~ year*log(impoth+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","impoth","year_impoth","imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:5]], 2, as.numeric)
agg <- aggregate(aggvars ~ hs10code, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code","year_mean","impoth_mean",
  "year_impoth_mean","imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:5] <- mm_df[,2:5] - mm_df[,7:10] 

mod1top10 <- lm(imp ~ year + year_impoth, data = mm_df)
mm_df$year_sector <- paste(newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2top10 <- lm(imp ~ year + year_impoth + year_sector, data= mm_df); 

## top 25
newdat <- data.frame(year = rep(c("20122016","20182020"), each = nrow(trdmat)), 
        hs10code = trdmat$hs10code, impoth = trdmat$imptop2520122016,
        imp = c(trdmat$impchn20122016/5, trdmat$impchn20212023/3)/1000
)

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_top25 <- newdat2
mod <- lm(log(imp+1) ~ year*log(impoth+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","impoth","year_impoth","imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:5]], 2, as.numeric)
agg <- aggregate(aggvars ~ hs10code, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code","year_mean","impoth_mean",
  "year_impoth_mean","imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:5] <- mm_df[,2:5] - mm_df[,7:10] 

mod1top25 <- lm(imp ~ year + year_impoth, data = mm_df)
mm_df$year_sector <- paste(newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2top25 <- lm(imp ~ year + year_impoth + year_sector, data= mm_df); 

## ALL
newdat <- data.frame(year = rep(c("20122016","20182020"), each = nrow(trdmat)), 
        hs10code = trdmat$hs10code, impoth = trdmat$impALL20122016,
        imp = c(trdmat$impchn20122016/5, trdmat$impchn20212023/3)/1000
)

postrade <- trdmat$hs10code[trdmat$impchn20122016/5 > 0 | trdmat$impwor20122016/5 > 0]
newdat2 <- newdat[newdat$hs10code %in% postrade,]
newdat2_ALL <- newdat2
mod <- lm(log(imp+1) ~ year*log(impoth+1), data = newdat2)
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs10code = newdat2$hs10code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","impoth","year_impoth","imp","hs10code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:5]], 2, as.numeric)
agg <- aggregate(aggvars ~ hs10code, FUN = mean, data = mm_df)
colnames(agg) <- c("hs10code","year_mean","impoth_mean",
  "year_impoth_mean","imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:5] <- mm_df[,2:5] - mm_df[,7:10] 

mod1ALL <- lm(imp ~ year + year_impoth, data = mm_df)
mm_df$year_sector <- paste(newdat2$year, substr(newdat2$hs10code,1,2),sep = "_")
mod2ALL <- lm(imp ~ year + year_impoth + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_top10$hs10code, 1, 6)
mods_top10 <- list(coeftest(mod1top10, vcov = vcovCL, cluster = hs6codes)[1:3,1:4], 
  coeftest(mod2top10, vcov = vcovCL, cluster = hs6codes)[1:3,1:4])

vars <- c(rownames(mods_top10[[1]])[c(3)])
tab <- matrix(data = NA, ncol = length(mods_top10), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top10)){mod <- mods_top10[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$Alternative markets total trade") 
N <- c(sum(summary(mod1top10)$df[1:2]),sum(summary(mod2top10)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab2123_top10 <- tab
ptable(cbind(rownames(tab), tab1820_top10, tab2123_top10), "paper/tables/TableA1top.tex")

hs6codes <- substr(newdat2_top25$hs10code, 1, 6)
mods_top25 <- list(coeftest(mod1top25, vcov = vcovCL, cluster = hs6codes)[1:3,1:4], 
  coeftest(mod2top25, vcov = vcovCL, cluster = hs6codes)[1:3,1:4])

vars <- c(rownames(mods_top25[[1]])[c(3)])
tab <- matrix(data = NA, ncol = length(mods_top25), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top25)){mod <- mods_top25[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$Alternative markets total trade") 
N <- c(sum(summary(mod1top25)$df[1:2]),sum(summary(mod2top25)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab2123_top25 <- tab
ptable(cbind(rownames(tab), tab1820_top25, tab2123_top25), "paper/tables/TableA1mid.tex")

hs6codes <- substr(newdat2_ALL$hs10code, 1, 6)
mods_ALL <- list(coeftest(mod1ALL, vcov = vcovCL, cluster = hs6codes)[1:3,1:4], 
  coeftest(mod2ALL, vcov = vcovCL, cluster = hs6codes)[1:3,1:4])

vars <- c(rownames(mods_ALL[[1]])[c(3)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$Alternative markets total trade") 
N <- c(sum(summary(mod1ALL)$df[1:2]),sum(summary(mod2ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))
tab2123_ALL <- tab
ptable(cbind(rownames(tab), tab1820_ALL, tab2123_ALL), "paper/tables/TableA1bot.tex")

